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Abstract 



We present the first non- mean- field calculation of the magnetization M(T) 
of YBa2Cu307_,5 both above and below the flux-lattice melting temperature 
T m {H). The results are in good agreement with experiment as a function of 
transverse applied field H . The effects of fluctuations in both order parameter 
-0(r) and magnetic induction B are included in the Ginzburg-Landau free 
energy functional: ^(r) fluctuates within the lowest Landau level in each layer, 
while B fluctuates uniformly according to the appropriate Boltzmann factor. 
The second derivative (d 2 M/dT 2 )n is predicted to be negative throughout 
the vortex liquid state and positive in the solid state. The discontinuities in 
entropy and magnetization at melting are calculated to be ~ 0.034 ks per flux 
line per layer and ~ 0.0014 emu cm -3 at a field of 50 kOe. 
PACS numbers: 74.60.-w, 74.40.+k, 74.25.Ha, 74.25.Dw 
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Even more than perfect conductivity, the magnetization M is a sensitive probe of the 
superconducting state. Type-I superconductors have an abrupt jump in magnetic suscepti- 
bility from to — 1/Att at the superconducting transition. Low-T c type-II superconductors 
have a more gradual behavior: the magnetization sets in continuously at the upper critical 
field H c2 (T), reaching —H/4ir only at the lower critical field H c i(T). In high-T c supercon- 
ductors, M also turns on continuously, but deviates from mean-field predictions. M may 
also be affected by the flux-lattice melting transition which occurs well below the mean- 
field superconducting transition temperature. To predict these deviations from mean-field 
behavior is a stringent test of any theory of fluctuations in the superconducting state. 

In this paper we calculate the magnetization of the most studied high-T c material, 
YBa 2 Cu 3 7 _5, including the effects of fluctuations. To our knowledge, this is the first 
non-mean-field calculation of M in YBa2Cu30y_5. We consider an applied field H parallel 
to the c axis, and also calculate the first-order flux-lattice melting curve T m (H) for the same 
material. Both M and T m (H) are in very good agreement with experiment over a range of 
magnetic fields. 

Our approach is to start from a Ginzburg-Landau free energy functional which includes 
the field energy in the form [FJ 
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Here ip(r) is the complex order parameter, A(r) the vector potential, B = V x A, e* = 
2e is the charge of a Cooper pair, H = Hz is the applied magnetic field intensity, and 
a(T,z), b, and m* are phenomenological parameters. The periodic ^-dependence of a(T,z) 
is introduced to characterize the underlying layered structure. We assume that spatial 
modulation of a(T, z), which acts as a potential for Cooper pairs, results in a tight-binding 
form for ijj(r) along the z-axis, thus creating the effective mass anisotropy 7. The integral is 
to be carried out over a fixed sample volume V = (0o/ H)N^N z s, where s is the periodicity 
of the layered structure, 0o = hc/2e is the flux quantum, N$ is the number of flux lines 



threading the sample, and N z is the number of layers. The Gibbs free energy density 
appropriate to an experiment at fixed T and H is 

g = -(k B T/V) lnTr^A exp(-G[V, A)/k B T). (2) 

We also define free energy 'per flux line per layer' according to Q$ = (4> s/H)Q; internal 
energy and entropy are defined in the same fashion. 

To determine the magnetization, we assume a uniformly fluctuating magnetic induction 
V x A(r) = Bz, so that B can in principle assume different values in the flux liquid and 
solid phases (but remains parallel to the c axis). This approximation is most reasonable 
in the extreme type-II limit (k> 1) when density correlations among the vortices are well 
developed and vortices are clearly defined, as expected in the neighborhood of the melting 
line; it may deteriorate very near the upper critical field H c2 (T). In addition, the assumption 
of uniform B||c should be best in the vortex liquid phase, where contributions to the local 
field B(r) come from many positionally uncorrelated vortex line segments, as suggested by 
Brandt |3J], but may still be reasonable in the vortex solid phase, even though B in that 
phase should develop a nonzero variance a 2 = (B(r) 2 ) — (B(r)) 2 . 

We evaluate the statistical average (2) using the following procedure. First, we expand 
the order parameter ip(r) in a basis which consists of products of lowest Landau level states 
of the operator (2m*) -1 (— ihV± — e*A/c) 2 in the ab plane, and Wannier functions from the 
lowest band of states of the operator a(T, z) — (2m*)~~ 1 h 2 d 2 / dz 2 in the c direction ||. This 
leads to an explicit form for all but the last term of the integrand in eq. (1), in terms of the 
complex coefficients Ck >n of the expansion, corresponding to the k th lowest Landau state in 
the n th layer ||. In order to impose periodic boundary conditions in the ab plane, we use 
the Landau gauge 0. Then the statistical average implied by eq. (2) is carried out by a 
Monte Carlo (MC) procedure, in which the coefficients Cfc >n , and also the average magnetic 
induction B, are considered as fluctuating MC variables. The entire procedure is closely 
analogous to the "constant pressure" MC ensemble well known in classical fluids J7|. The 
variables analogous to pressure and volume are H and B. The MC step which changes B 
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increases or decreases the area of the ab plane at constant vortex number. 

Within this approach, the mean- field approximation to (2) is specified by a pair ip(r) 
and B, for which G is minimum ||. In the normal state (H > H^(T)) this is achieved for 
•0(r) =0 and B = H. In the mixed state, (H < H&^T)), the corresponding minimum is 
attained for a triangular lattice of straight vortex lines and magnetization 

B-H H-H c2 (T) 
~ 4vr 47r(2/t 2 /3 A - 1) ' K ) 

where (3a = 1.15959. . . is the Abrikosov ratio. In the limit k ^> 1 studied here, this for- 
mula becomes identical to the original Abrikosov result, which has 4tt(2k 2 — 1)(3a in the 
denominator. Thus our approach and approximations are indeed correct at the mean- 
field level in this limit. The corresponding mean- field free energy density is Q MF = 
-(87r)- 1 [H c2 (T)-H}y(2K'f3 A -l). 

In order to execute the MC calculation for YBa2Cu30y_5, we use the following set of 
parameters: T c0 = 93 K, dH c2 (T)/dT = -1.8 x 10 4 Oe/K, s = 11.4 A, 7 = 5, and k = 52. 
Although there is some experimental evidence that k varies with magnetic field |J, we 
neglect this field-dependence here. In most of our calculations, we have considered a cell 
containing = 10 2 vortices and N z = 10 layers. Our results are based typically on 2-3 xlO 5 
MC passes through the entire sample following ~2x 10 4 MC passes for equilibration. 

Fig. 1 shows the calculated magnetization M(T) of YBa2Cu30y_5 as a function of tem- 
perature T at three different values of H. The mean-field predictions are shown for compar- 
ison. There is no sign of a true phase transition at the mean-field transition temperature 
T c2 (H), since (dM/dT) H is continuous at that point. Instead, there is an apparently first- 
order melting transition at a lower temperature T m (H), as signaled by a weak discontinuity 
in the magnetization curves (denoted by arrows in the Figure). The melting curve T m (H) 
inferred from this discontinuity is shown in the inset to Fig. 1; it is in good agreement with 
experiment. 



The general behavior of M(H, T) agrees very well with experiment [IT]. For example, the 



calculated second derivative (d 2 M/dT 2 )H is negative throughout the vortex liquid phase, in 



accord with recent measurements based on a differential torque technique ||12|| . A second- 
degree polynomial fit to our magnetization results just above the melting temperature T m 
yields (d 2 M/dT 2 ) H = (-0.0038±0.002) emu cm^K" 2 for H = 20 kOe, and (d 2 M/dT 2 ) H = 



(—0.0030 ± 0.002) emu cm _3 K~ 2 for H = 50 kOe. Experiment jTIJ also gives a negative 
(d 2 M/dT 2 )H for fields in the range 10-20 kOe, and of the same order of magnitude. In the 
solid phase, we find our calculated (d 2 M/dT 2 )H > 0. 

Another striking feature of our results is the crossing of the magnetization curves as 
a function of temperature. This crossing is also observed in experiment at a similar tem- 
perature ||11||. In E^S^CaC^Os+a, the same phenomenon has attracted much theoretical 



attention |T3|Jl^| , because it is believed to occur at a unique temperature T* independent 
of field. In YBa2Cu307_,5, experimental data of Welp et al. [|llj reveal that magnetization 
curves do not cross at a single point, presumably because, although layered, YBa 2 Cu 3 7 _5 
is only moderately anisotropic. As may be seen in the Figure, our calculated magnetization 
curves also fail to cross at a single point. Furthermore, the order of the calculated crossings 



as a function of field is the same as observed in YE^CusCv-^ [IT|. 



The second derivative of the magnetization can also be calculated using a Maxwell rela- 
tion 

where Ch = —T[d 2 Q jdT 2 )u is the specific heat of the sample at constant H |L5[]. Since 
eq. (4) equates the second derivative of one macroscopic quantity to the first deriva- 
tive of another, it provides a potentially more precise way of determining (<9 2 M/<9T 2 )# 
than a direct calculation of M(T). Note that in the mean-field approximation Cf} F = 
-T(d 2 G MF /dT 2 ) H = (T/4n)(dH c2 /dT) 2 /(2K 2 f3 A - 1), which is independent of H. This im- 
plies, via eq. (4), that M MF (T) is linear in T, in agreement with eq. (3). To go beyond mean- 
field theory, we have done an MC calculation of Ch{H,T), using the fluctuation-dissipation 
theorem |TB[. Our results are presented in Fig. 1 (inset). Clearly, d(Cn/T) / dH < in 
the vortex liquid phase, which implies, in agreement with experiment fl2|], that M(T) has 



negative curvature throughout the liquid phase, not just in the vicinity of the mean-field 
critical temperature T c2 (H). The straight lines in the Figure are constructed using slopes 
determined from independent calculations of the magnetization. To an excellent approxi- 
mation, they are tangent to the specific heat curves, as expected from eq. (4), thus con- 
firming the self-consistency of our approach. Experimental specific heat data |T7] confirm 



that oI(Ch /T)/dH < in the vortex liquid phase. In the solid phase, our work predicts 
that d(CH/T)/dH > 0. This is in agreement with recent, high-accuracy specific heat mea- 
surements by Schilling and Jeandupeux [|l^] on large twinned YBa2Cu307_,5 crystals. The 
interesting region in the immediate neighborhood of the melting line has not been resolved, 
however. 

Fig. 2 shows the in-plane structure factor S^qj^O) defined as the thermal average 

5(q±,0) = dWr^(r)|V(0|V q - (r - r,) ), (5) 

for a field H = 50 kOe at two temperatures T = 82.8 K and T = 83.0 K, corresponding 
to the vortex solid and vortex liquid phases. This dramatic change occurs at the temper- 
ature of the magnetization discontinuity, thus confirming that this discontinuity signals a 
melting transition. In each case, the displayed S^qj^O) is the result of averaging over 100 
configurations from the thermal ensemble. For clarity, we have removed the central maxi- 
mum at q_|_ = 0, and have divided by the "atomic" structure factor exp(— qj_£ 2 /4), where 
i = (0O/27T5) 1 / 2 is the magnetic length. The regular periodic structure of the maxima in 
Fig. 2(a) corresponds to an ordered crystalline phase, while the concentric rings of Fig. 2(b) 
characterize an isotropic fluid. In both cases, the deviations from perfect isotropy can be 
attributed to the finite size of the sample, as well as to its rectangular shape. 

Finally, we address the issue of the order of the melting transition. As indicated by Fig. 2, 
the vortex solid state is a phase of discrete symmetry, whereas the vortex liquid is an isotropic 
phase of much higher symmetry. Therefore, upon melting, the vortex ensemble undergoes 
a discontinuous symmetry change. By the well-known argument of Landau [TJJ, the vortex 
melting transition has to be first-order, with a finite jump in entropy per unit volume AS. 
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A critical point (defined by AS = 0) cannot exist, and the liquid-solid phase boundary must 
terminate by intersecting either the coordinate axes or other phase boundaries. 

To calculate AS, we use a variant of the histogram method of Lee and Kosterlitz |2U 



In principle, one should resolve the energy distribution P(G) into two Gaussian peaks at 
T m , then confirm that the dip between the peaks scales like a surface energy with increasing 
sample size. Such a calculation was done in Ref. for the 3D XY model. In the present 



case, this would be a formidable computational effort, because AS is small and the two 
peaks are not resolved at any accessible system sizes. Instead, we perform a long MC run 
(~ 10 6 passes through the entire lattice) near T m , starting from the system ground state. 
During the simulation, the system flips ~ 2-4 times between the two states in equilibrium. 
In our standard diffusive sampling algorithm, we can identify continuous (in MC "time") 
sequences of representatives belonging to the same homogeneous phase. The two phases can 
be clearly distinguished by their structure factors and mean internal energy. 

In Fig. 3 we plot the probability distribution of internal energy, P(G), for these two 
phases in equilibrium at T m ~ 83 K and H = 50 kOe, at three system sizes. The low-energy 
peak always corresponds to the ordered vortex solid phase. Since the entropy change AS^ 
decreases with system size, the value AS^ ~ 0.034 ks should be considered an upper bound 
to AS^ in the thermodynamic limit. AS^H = 20 kOe) has a similar size dependence and 
is ~ 30% smaller than AS^(H = 50 kOe). This is the expected trend, since T m is higher 
at H = 20 kOe, and typical superconducting parameters, such as the condensation energy 
G MF , decrease with increasing temperature. 

From our calculated AS at melting and the computed slope (dH/dT) m of the melting 
curve from Fig. 1, we can estimate the magnetization jump AM at melting via the Clausius- 
Clapeyron relation AS/ AM = —(dH/dT) m . Inserting the calculated values of this slope and 
of AS^, we obtain AM ~ 0.0014 emu cm" 3 at H = 50 kOe, and AM ~ 0.0005 emu cm" 3 at 
H = 20 kOe. These values are consistent with the directly calculated AM seen in Fig. 1. In 
experiment, no finite magnetization jump has yet been observed |T2"| , |2"2"|] , although transport 
measurements on untwinned YB^CusOts crystals are widely interpreted as evidence for a 



first-order melting transition |10|j23[| . The null result of Farrell et al. j[2] for AM puts an 
upper bound AS^ < 0.003 ks at H = 20 kOe, a value almost ten times smaller than predicted 
in the present work. However, as suggested by these authors themselves, the presence of 
defects may have suppressed the measured AM to some extent. This discrepancy remains 
to be settled by an experiment which attains perfect reversibility in both temperature and 
field. 

To summarize, we have presented the first non-mean-field calculation of magnetization 
in YBa2Cu307_,5, using a constant-H Monte Carlo technique in conjunction with a lowest 
Landau level approximation. Our results yield a magnetization in very good agreement 
with experiment in both the flux liquid and flux lattice state, as well as a melting curve 
very close to experiment. Our results provide perhaps the most detailed evidence to date 
that a Ginzburg-Landau free energy functional, based on a complex scalar order parameter, 
adequately describes the thermodynamic properties of YB^CusOy-^ near the flux lattice 
melting curve. 

We are grateful for valuable conversations with Professors D. E. Farrell, T. R. Lemberger, 
and W. F. Saam. This work was supported by NSF through Grant DMR94-02131, the 
Midwest Superconductivity Consortium through DOE Grant DE-FG02-90ER45427. The 
calculations were performed on the local network of DEC Alpha stations. 
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FIGURE CAPTIONS 



1. Calculated magnetization M(T) of YBa 2 Cu 3 7 _ ( 5 at three different applied fields, 
H = 10 kOe, 20 kOe, and 50 kOe. Dashed lines represent the mean- field solution 
(3). Solid lines are spline curves connecting the calculated points. Arrows denote 
melting temperatures T m (H), as determined by the calculated discontinuity in M(T). 
Estimated errors in M(T) are much smaller than the symbol sizes. x N z = 1 2 x 10. 
Right inset: locus of the liquid-solid phase boundary in the H- T plane, as determined 



by our calculations and as measured by Farrell et al. |§ and Safar et al. [10 |. 
Left inset: specific heat Ch as a function of magnetic field, taken at two temperatures, 
T = 83 K and 87 K. The dashed line represents the mean-field value C^ F /T. Arrows 
indicate the approximate location of the fields at melting H m (T). Straight lines have 
slopes -0.0038 emu cm" 3 K~ 2 at 87 K and -0.0030 emu cm" 3 Kr 2 at 83 K (see text), 
x N z = 6 2 x 6. 

2. (a) In-plane structure factor S(qj_,0), divided by the "atomic" structure factor 
exp(— qj_£ 2 /4), taken in the solid phase at T = 82.8 K and H = 50 kOe, and av- 
eraged over 100 configurations. The central maximum at = has been removed 
for clarity. 

(b) Same as (a), but in the liquid phase at T = 83.0 K. 

3. Probability distribution P(G) of internal energy at the finite-size melting point T m ~ 
83 K, H = 50 kOe, and three different system sizes. For each system size, the two 
separate peaks represent the distributions of energy within the solid and the liquid at 
the same T and H. Because of the weak dependence of T m on system size, we have 
used for the zero of a size-dependent energy. Horizontal bars denote the calculated 
latent heats per flux line per layer, in units of ksT. 
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